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INTRODUCTION 


Time  Delay  Estimation  (TDE)  has  received  considerable 
attention  in  recent  years,  primarily  because  of  its 
applicability  to  source  position  estimation  in  passive 
sonar.  This  paper  presents  the  more  significant  related 
theory  and  analysis  developed  recently  in  this  area.  The 
two  broad  classes  of  (1)  TDE  when  the  source  and  receivers 
are  in  relatively  fixed  positions  and  (2)  TDE  when  there  is 
relative  motion  between  the  source  and  receivers  are 
considered,  with  emphasis  on  the  latter,  more  general  class. 
By  way  of  background,  a  brief  description  of  the  typical 
sonar  problem  which  has  motivated  much  of  the  research  in 
this  area  is  given.  Next  the  no  motion  case  is  studied, 
followed  by  analysis  for  the  case  with  motion.  Finally,  a 
computer  implementation  of  a  motion  compensated 
cross-correlator  is  presented  and  results  which  corroborate 
the  theory  and  analysis  are  presented  and  discussed.  A 
brief  computer  program  description  is  provided  in  Appendix 
A. 


1.0  APPLICATION  OF  TIME  DELAY  ESTIMATION  TO  SONAR 


Sonar  systems  may  be  divided  into  two  major  categories: 
Active  Sonar  and  Passive  Sonar.  In  Active  Sonar,  a 
pre-determined  signal  is  transmitted  by  the  sonar  into  the 
ocean  environment.  The  received  echo  from  objects  in  the 
ocean  is  then  analyzed  to  obtain  information  about  the 
objects  such  as  their  nature,  location,  and  movement.  Since 
the  form  of  the  signal  can  be  controlled,  it  can  be  designed 
to  best  meet  some  desired  performance  criteria  such  as 
detectability  and  spatial  resolution.  The  distance  to  an 
object  is  determined  by  measuring  the  total  travel  time  of 
the  signal,  from  the  transmitter  to  the  object  and  back  to 
the  sonar  receiver,  and  utilizing  knowledge  of  the  speed  of 
sound  in  water.  The  angle  to  the  object  can  be  determined 
by  processing  the  received  data  only  in  certain  angular 
directions  (a  process  referred  to  as  beamforming).  The 
major  drawback  of  Active  Sonar  is  that  the  location  of  the 
sonar  is  given  away  when  it  transmits,  a  very  undesirable 
effect  for  military  applications.  Also,  its  performance  can 
be  greatly  limited  by  reverberation  of  the  transmitted 
signal,  especially  for  high-powered  transmitters  in  shallow 
water  or  in  an  area  with  many  scatterers,  such  as 
particulate  matter  and  biologies. 

The  alternative  to  Active  Sonar  is  Passive  Sonar  which 
"listens"  for  sounds  emanating  from  the  objects  themselves, 
such  as  flow  noise,  machinery  noise,  and  marine  life  noises. 
The  signal  is  typically  processed  by  spectral  analysis  to 
determine  the  characteristics  of  the  objects.  A  major 
limitation  to  Passive  Sonar  is  a  lack  of  complete  knowledge 
of  the  characteristics  of  the  signal,  requiring  more  flexi¬ 
bility  in  the  processing.  Additionally,  since  the  travel 
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time  of  the  received  signal  is  not  known,  the  distance  to 
the  source  cannot  be  measured  directly.  Instead,  multiple, 
accurate  bearings  to  the  source  must  be  measured,  from  which 
a  range  can  be  determined  indirectly.  It  is  primarily  for 
the  purpose  of  obtaining  these  accurate  bearings  that  Time 
Delay  Estimation  (TDE)  is  employed  in  Passive  Sonar. 
References  [1],  [2],  [3],  and  [4]  contain  extensive  material 

on  sonar  and  related  signal  processing.  Reference  [5] 
contains  many  papers  on  TDE  and  provides  an  excellent 
overview  of  much  of  the  more  recent  work  on  the  subject;  in 
particular,  the  article  by  Carter  summarizes  TDE  as  applied 
to  Passive  Sonar. 

Figure  1  illustrates  the  basic  model  used  for  analysis 
of  passive  TDE.  The  source  and  two  receivers,  with  a  single 
propagation  path  to  each  receiver,  are  all  assumed  to  lie  in 
the  same  plane.  In  this  planar  model  it  is  also  assumed 
that  there  is  no  relative  motion  among  the  source  and 
receivers.  The  signal,  s(t),  arrives  at  Receiver  1  then  at 
Receiver  2  at  a  time  D  later.  For  a  distant  source  the 
attenuation,  <*,  of  the  signal  arriving  at  Receiver  2  is  =  1. 
Also  input  at  each  receiver  are  noise  n^t)  and  n2(t),  which 
are  assumed  to  be  uncorrelated  with  each  other  and  the 
signal.  The  signal  and  noise  are  also  assumed  to  be 
zero-mean,  wide-sense  stationary,  broadband  random 
processes.  The  equations  for  the  model  are 

rx(t)  =  s  ( t )  +  n^t)  (1.1. a) 

r2  ( t)  =  s(t  +  D)  +  n2(t)  (1.1. b) 

There  are  more  than  one  possible  source  locations  for 
which  the  same  time  delay  value,  D,  would  occur  between 
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Sou  r  c  e 


FIGURE  1 
PLANAR  MODEL 
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arrivals  at  the  two  receivers  [6],  as  illustrated  in  Figure 
2.  In  fact,  the  source  could  be  anywhere  on  the  locus  of 
points  that  form  a  hyperbola.  Figure  3  shows  the  geometry 
from  which  the  equation  for  the  hyperbola  and  the  equation 
for  the  bearing  to  the  source  (measured  from  the  midpoint 
between  the  receivers)  for  any  point  on  the  hyperbola  can  be 
derived.  The  difference  in  the  range  between  the  source  and 
Receiver  1  and  Receiver  2  is 


R2  -  R^  =  cD 


where  c  is  the  speed  of  sound  in  water,  D  is  the  time  delay. 

Substituting  for  R^  and  R, ,  the  equivalent  Euclidian 
distance  for  the  given  coordinates 


•yj (X+L/2  )2  +  Y 
(4/c2D2)X2  - 


L/2 )' 


[4/  (4(L/2  )" 


+  Y2  =  CD 
-  c2D2 ) ] Y2 


1. 


(1.2) 


Equation  (1.2)  defines  a  hyperbola.  The  bearing  angle  0  is 
found  from 

COS ( 0 )  =  ( cD/L ) [ 1  +  (L/2R)2  -  ( CD/2R)2 ] 1/2 

As  seen  in  Figure  4,  if  the  distance  to  the  source  is  large 
compared  to  the  receiver  separation  (i.e.,  R>>L,  as  is 
usually  the  case),  then  the  bearing  is  approximated  by  the 
angles  formed  by  two  straight  lines  which  form  assymptotes 
to  the  hyperbola.  This  bearing  then  is  related  to  the  time 
delay  D,  and  the  sensor  separation,  L,  by 
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FIGURE  2 

MULTIPLE  SOURCE  POSITIONS  WITH  SAME  DELAY 
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FIGURE  3 

SOURCE/RECEIVER  GEOMETRY 
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FIGURE  4 

BEARING  ANGLE  FOR  LARGE  RANGE 
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COS(0)  3  (CD/L) 

0  £  COS-1  (CD/L)  (1.3) 

where  c  is  the  speed  of  sound  in  the  ocean.  Thus,  for  known 
sensor  spacing,  L,  an  estimate  of  the  time  delay  yields  the 
bearing  to  the  source  (albeit,  it  is  an  ambiguous  bearing 
which  must  be  resolved  by  other  means  such  as  ship 
maneuvers,  or  auxilliary  information  like  rough  bearings 
from  passive  beamforming).  If  three  co-linear  receivers  are 
used,  then  ranging  is  possible  by  using  two  bearings  and  the 
sensor  spacing  for  a  triangulation  [5,  pp.  463-470]. 


2.0  TIME  DELAY  ESTIMATION  WITH  FIXED  SOURCE  AND  RECEIVERS 


The  most  common  method  for  estimating  time  delay 
passively  is  to  compute  the  cross-correlation  of  the  signal 
from  two  spatially  separated  receivers.  A  generalized 
method  for  estimating  time-delay  with  cross-correlation  was 
described  by  Knapp  and  Carter  [7],  and  is  recounted  here. 


For  the  model  defined  by  (1.1) 


Rrlr2(tl‘ 


* 


[rl(tl)r2(t2)] 


The  value  of  x  =  -  tj  which  gives  the  maximum  value  for 

the  cross-correlation  is  chosen  as  the  estimate  of  the  true 
delay,  D. 


Since  we  are  always  limited  to  a  finite  observation 
time,  Rrir2(T)  must  be  estimated.  For  ergodic  processes, 
this  is  computed  by: 


R  _ 
rlr2 


(i) 


~  r1(t)r2(t-t )dt 


(2.1) 


where  T  is  the  finite  observation  time.  This  estimation  of 

A 

Rrlr2^T^  P^us  the  presence  of  noise  will  result  in  errors  in 
the  computation  of  D,  the  delay  estimate.  In  practice, 
Rrir2^t)  is  computed  from  the  inverse  Fourier  transform  of 
the  estimate  of  the  cross-spectrum,  Grir2^^'  by: 

Rrlr2(T)  =  sZ  Grlr2(f)  ej2nfT  df  (2.2) 

In  order  to  reduce  the  error  in  Rrir2(T)'  an<^  bbus  improve 
the  accuracy  of  D,  we  can  pre-filter  r^(t)  and  r^(t)  prior 
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to  cross-correlating.  The  pre-filters  would  be  selected  to 
emphasize  frequencies  where  the  signal-to-noise-ratio  (SNR) 
is  highest  and  de-emphaize  them  where  it  is  low.  Figure  5 
illustrates  such  a  configuration. 

For  pre-filters  H^(f)  and  H2(f),  which  operate  on 
r1(t)  and  r2(t)  and  give  outputs  y1(t)  and  y2(t), 
respectively,  the  cross-spectrum  is  then 

Gyly2  (f)  =  Hl(f>  H2(f)*  Grlr2(f)  (2.3) 

and  the  cross-correlation  is  computed  from 

Ryly2  (t)  =  C  w(f>  Grlr2(f)  df  (2-4) 

where  W(f)  =  H^(f)  H2(f)  .  Thus,  the  prefiltering  can  be 
done  as  frequency  weighting  in  the  frequency  domain. 

Equation  (2.4)  is  referred  to  as  the  Generalized 
Cross-Correlation  (GCC)  function,  and  W(f),  the  General 
Frequency  Weighting,  which  is  chosen  to  optimize  some 
performance  criterion. 

2 • 1  Some  Frequency  Weightings  Used  In  Cross-Correlation 

Table  1  lists  several  proposed  frequency  weightings. 
Each  is  described  in  detail  in  [7].  As  the  weightings  are 
functions  of  the  signal  and  noise  spectra,  it  must  be  kept 
in  mind  that,  in  practice,  the  weightings  themselves  must  be 
estimated  (unless  there  is  sufficient  prior  knowledge  of  the 
spectra).  The  major  features  of  each  are  outlined  in  the 
following. 
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FIGURE  5 

GENERALIZED  CROSS-CORRELATION  FUNCTION 
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PROCESSOR 


CROSS-CORRELATION 


ROTH 


SCOT 


PHAT 


ECKART 


ML  or  HT 


WEIGTH:  W(f)  =  H1(f)H2*(f) 


1 

1 


Grlrl 


(f) 


1 


■yGrlrl(f)  Gr2r2(f* 


1 


Grlr2(f ^ 


tGnlnl<f>Gn2n2<f>J 

Crlr2<f> 


Grlr2(f) 


t1-cnr2<f)] 


TABLE  1 

SOME  FREQUENCY  WEIGHTING  FUNCTIONS  FOR  CROSS-CORRELATION 
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The  first  weighting,  W(f)  =  1  is,  of  course,  no 
weighting  at  all  and  is  used  when  the  signal  and  noise  have 
flat  spectra,  equal  bandwidths  and  the  SNR  is  constant 
across  the  entire  band.  It  is  often  referred  to  as  the 
Standard  Cross-Correlation  (SCC). 

The  Roth  processor  [8],  which  gives  correlation 
estimate 


Grlr2(f J 


RylY2(t)  =  L 


2nft 

e  df 


Grlrl<f> 


is  also  the  estimate  of  the  impulse  response  of  the  filter 


Grlr2(f) 
H(f)  =  - 

Grlrl<f> 


which  is  an  optimum  approximation  of  the  linear  mapping  of 
r1(t)  to  r2(t)  (Wiener-Hopf  filter).  Since 


Grlrl<f>  =  Gss<f> 


Gnlnl(f) 


we  get 


Ryly2 ^ x  ^  I® 


Grlr2(f ^ 


27tfT 

e  df 


Gss<f>  +  Gmnl'f> 


This  weighting  can  be  seen  to  suppress  those  frequencies 
where  the  noise  at  receiver  1  is  highest. 
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The  Smooth  Coherence  Transform  (SCOT)  [9],  [10] 

processor  can  be  seen  to  suppress  frequencies  where  the 
noise  in  either  or  both  channels  is  high.  The  process  is 
equivalent  to  "pre-whitening"  the  inputs  r^(t)  and  r2 (t)  by 

Hx(f)  =  VGrlrl(f)  and  H2(f)  =  1/Gr2r2(f) 

and  then  cross-correlating  the  results.  For  the  case  of 
white  signal  plus  non-white  noise,  the  effect  of 
pre-whitening  is  to  attenuate  the  cross-spectrum  estimate 
where  the  noise  is  highest  and  conversely  where  it  is 
lowest.  If  Gr-^r^(f)  =  Gr2r2^^'  SC0T  can  seen  to 

equivalent  to  the  Roth  processor. 


The  Phase  Transform  (PHAT)  [10]  window  was  developed 
specifically  to  minimize  the  spreading  of  the  correlation 
peak  caused  by  non-white  signal.  For  the  ideal  case  of 
white  signal  and  no  noise,  the  cross-correlation  is  a  delta 
function,  6(t  -  D).  A  non-white  signal  results  in  the 
spreading  or  smearing  of  the  delta  function  through 
convolution  with  the  transform  of  the  signal  spectrum.  For 

-A 

perfect  estimation,  where  Grlr2(f)  =  Grlr2(f),  then 


Grlr2<f> 


Prlr2^f ) 


0j2nfD  _  ej6 (f ) _ 


That  is  we  get  unit  magnitude  with  linear  phase.  Then 
Ryly2(t)  =  5 (t  -  D)  and  no  spreading  takes  place.  In 
practice,  though,  Grlr2(f)  ?  Grlr2(f)'  and  e(f)  ?  2nfD,  but 
is  erratic,  resulting  in  a  correlation  which  is  not  a  delta 
function.  Another  serious  problem  with  the  PHAT  is  that  it 
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weights  Grlr2(f)  as  the  inverse  of  G  (f),  thus 
accentuating,  rather  than  attenuating  the  frequency  regions 
where  the  signal  is  weakest. 


The  weighting  referred  to  as  the  Eckart  Filter  [11]  has 

the  characteristic  of  maximizing  the  ratio  of  the  change  in 

the  mean  of  the  correlator  output  when  signal  is  present,  to 

the  standard  deviation  of  the  correlator  output  when  there's 

noise  alone.  Like  the  SCOT,  this  processor  suppresses 

frequencies  where  the  noise  is  highest,  but,  additionally, 

it  gives  zero  weighting  to  the  regions  where  G  (f)  =  0. 

s  s 

The  last  weighting  function,  the  HT  (for  the 
originators  Hannan  and  Thompson  [12])  is  shown  in  the  next 
section  to  be  a  maximum  likelihood  estimator  of  the  delay, 

D,  when  the  assumption  of  Gaussian  noise  and  signal  is  made. 
It  is  also  shown  that  the  variance  of  the  delay  estimate, 
Var[D]  is  equal  to  the  Cramer-Rao  Lower  Bound  ( CRLB ) ,  and, 
as  such,  it  is  an  optimum  (minimum  variance)  estimator  (see 
next  section).  It's  use  of  the  magnitude  squared  coherence, 


Grlr2<f> 


Crlr2(f)  G 


rlrl(f)  Gr2r2^f) 


(2.5) 


between  the  receivers  serves  to  emphasize  the  high  SNR 
regions  while  attenuating  the  low,  like  the  SCOT.  Like  the 
PHAT,  the  Grir2^^|  term  in  the  denominator  tends  to 
reduce  the  spreading  of  the  correlation  peak,  but  in  this 
case,  it  suppresses  the  regions  where  the  phase,  6(f),  is 
most  erratic,  (namely  the  low  or  zero  SNR  regions)  by  the 
coherence  terms . 
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For  the  case  where  the  SNR  is  low  and  «  =  1,  the  HT  and 
Eckart  processors  can  be  shown  to  be  equivalent  [7]. 
Additionally,  at  low  SNR,  they  could  be  interpreted  as  SCOT 
processors  with  added  SNR  weighting  or  as  PHAT  processors 
with  added  SNR  weighting.  These  processors  have  recently 
been  implemented  and  analyzed  as  described  in  [13],  [14], 

[15]  and  [16]. 

2.2  Derivation  of  Generalized  Cross-Correlation  (GCC) 

as  a  Maximum  Likelihood  Estimate  r— ^ 

Knapp  and  Carter  [7]  derive  the  Maximum  Likelihood  (ML) 
estimate  for  time  delay,  which  proves  to  be  the  same  as  that 
proposed  by  Hannan  and  Thompson.  This  derivation  is 
summarized  by  Scarbrough  [13]  with  some  added  clarification 
and  is  presented  here  in  a  slightly  modified  form. 

Outlining  briefly,  the  assumption  that  the  signal  and 
noise  are  Gaussian  processes  is  made  in  addition  to  the 
previous  assumption  of  being  wideband,  wide  sense 
stationary,  zero-mean,  and  uncorrelated  with  each  other.  We 
will  determine  an  observation  vector,  R,  which  depends  upon 
the  true  time  delay,  D,  and  the  signal  and  noise  auto-  and 
cross-power  spectra,  which  are  assumed  known.  A  maximum 
likelihood  (ML)  estimate  of  the  delay  is  made  by  picking  as 
our  estimate  the  hypothesized  delay,  t ,  which  maximizes  the 
probability  density  function  (pdf)  of  R  conditioned  on  t . 

It  will  be  seen  that  the  estimator  reduces  to  picking  the  t 
which  maximizes  the  cross-correlation  function  of  the 
receiver  signals  whose  cross-spectrum  has  been  multiplied  by 
a  weighting  function  formed  from  the  auto-  and  cross-spectra 
of  the  signal  and  noise.  Furthermore,  it  will  be  seen  that 
the  variance  of  this  ML  estimate  of  delay  achieves  the 


18 


Cramer-Rao  Lower  Bound  (CRLB)  and,  as  such,  is  the  optimum 
estimator,  in  a  minimum  variance  sense. 

Our  first  step  is  to  express  the  received  waveforms, 
r^(t)  and  ^(t)  over  the  interval,  0-»T,  by  their  Fourier 
Series  coefficients.  This  can  be  done  by  forming  the 
periodic  extension: 

ri(t+mT)  =  Ri(t),  O^tST,  i  =  1,  2 

for  any  integer  value  of  m.  The  Fourier  Series  coefficents 
for  these  periodic  extensions  are  then 

Ri(k)  =  1/T  Si  ri(t)e":ik27xfot  dt,  i  =  1,2  (2.6) 


where 


f0  =  l/T 

The  two  received  signals  r^(t),  i  =  1,  2  can  be  computed 
from  these  Fourier  coefficients  by  taking  one  period  of  the 
inverse  relationship.  That  is 
00 

?i( t)  =  ^Ri(k)e;ik27xfot,  for  all  t 
k=-» 

r^(t)  =  r^(t)  for  O^t^T 

For  practical  cases  the  r^(t)  are  bandlimited  signals  so 
that  the  Fourier  coefficients  are  all  zero  outside  some 
range,  say,  -N^kSN.  We  see  then  that  the  r- (t)  can  be 
characterized  by  the  set  of  Fourier  Series  coefficients, 
Ri(k),  -NSskSN. 
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We  wish  to  obtain  an  expression  for  the  joint 
probability  density  function  of  these  coefficients 
conditioned  on  a  hypothesized  value  for  the  time  delay, 
denoted  by  t,  which  is  chosen  from  the  range  of  possible 
delays.  As  will  be  seen,  the  coefficients  are  also 
dependent  upon  the  power  spectral  density  functions  of 
received  signals,  which,  for  now,  we  will  assume  are  known. 

First  we  observe  from  (2.6)  that  each  coefficient 
Rj_(k) ,  i  =  1,  2,  for  any  particular  value  of  k  in  the  range 
-NSkSN,  is  a  zero  mean  Gaussian  random  variable.  This  is 
because  it  is  derived  by  a  linear  transformation  on  the 
received  signal  r^(f),  which,  as  we  have  assumed,  is 
zero-mean  Gaussian.  So  we  have 

E[R1(k)]  =  0  -NSkSN 

E[R2(k)]  =  0  -NS kSN 

For  large  T  (i.e.,  small  f  ) 

o 


E[R1(k)R2(j)*]  = 


for  k  =  j 


for  k  f  j 


That  is  for  large  T,  the  coefficients  are  essentially 
uncorrelated. 

We  next  form  a  bivariate  Gaussian  random  column  vector 


R(k)  4  [R1(k),R2(k)] 
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for  each  value  of  k  in  -N^k^N  (  the  '  denotes  the  transpose 
of  the  vector).  In  order  to  represent  its  conditional 
density  function,  we  need  the  covariance  matrix 

Cov[R1(k) ,R1(k) ]  Cov[R1(k) , R2(k) ] 

Cov [R2 ( k ) , rx ( k ) ]  Cov [R2 ( k ) , R2 ( k ) ] 

where  Cov  (Rm/Rn)  is  the  covariance  of  the  random  variables 
Rm(k)  and  Rn(k) 

Cov[Rm,Rn]  =  E[Rm(k)Rn(k)*|T]  -  E[Rm(k)|t]E[Rm(j)*|t] 
Since  the  R^(k)  are  zero  mean  random  variables 
Cov[Rm,Rn]  =  E[Rm(k)Rn(k)*|t] 
which  from  (2.6)  is 

cov[Vk)  ,  Rn(k)]  =  |  <5rm,rn<kf0) 

Substituting  into  the  covariance  matrix: 


COV 


k|  t 


COV 


k  1 1  T 


Griri(kfo>  Grlr2(kfo) 


Gr2rl(kfo)  Gr2r2(k£o) 


=  f  Q(k) , 


(2.7) 


where  Q(k)  is  referred  to  as  the  cross-spectral  density 
matrix.  For  any  particular  value  of  k,  -N^k^N,  we  now  can 
write  the  bivariate  Gaussian  probability  density  function 
for  R(k)  conditioned  on  t ,  in  standard  form: 
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p(R(k)|x)  =  \2n 


cov. 


k|x 


1/2|  exp  -1/2  | 


R(k)*’ [COVklx]_1  R(k 


! 


We  now  wish  to  consider  the  conditional  pdf  of  R(k)  for  all 
values  of  k,  -N^k^N,  considered  jointly.  To  do  this  we 
define  our  observation  vector  as 


R=  [R ( -N) ,  R( -N+l ) ,  .  .  .  R(N) ] . 


Now  since  the  R(k)'s  are  uncorrelated  (and  being  Gaussian, 
also  independent),  the  conditional  joint  pdf  of  R  is  simply 
the  product  of  the  pdf's  for  each  individual  R(k): 

N 

P(RIt)  =  7T  p(R(k)  |  t  ) 

k=-N 

if  we  let: 


c  =  TT  (27t|c0Vk  |1/2) 

k=-N  ' 


where 


is  the  determinant,  and 


N 


J1  =  Z  R(k)*' [COVk)x]  XR(k) 


k=-N 


(2.8) 


(2.9) 


we  can  write  p  (R|x)  in  the  simple  form  of 


p(R| x )  =  c  exp(-J1/2 ) 


(2.10) 
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Summarizing  where  we  stand,  we  have  an  observation  vector  R 
which  is  conditioned  upon  the  possible  values  of  the  time 
delay,  and  is  also  dependent  upon  the  power-spectra  of  the 
received  data.  The  elements  of  our  observation  vector  are 
the  Fourier  Series  coefficients  of  r1(t)  and  r2(t).  Now  we 
will  use  a  maximum  a  posteriori  decision  criterion  and 
assume  that  all  possible  values  of  x  are  equally  likely  (the 
actual  range  of  possible  values  for  x  is  dependent  upon  the 
source/receiver  geometry).  Thus  the  maximum  likelihood 
estimate  for  D  (the  actual  time  delay  of  interest),  is  to 
choose  D  =  x ,  where  x  gives  the  maximum  value  for  p  (R  |  x ) . 
That  is,  we  choose  the  x  that  most  likely  would  have  given 
us  the  observed  values  in  R. 

We  will  now  use  equations  (2.7),  (2.8),  (2.9),  and 
(2.10)  and  some  simplifying  approximations  to  reduce  this 
test  to  equivalent  ML  tests  until  we  get  it  in  the  final 
desired  form. 


First,  we  show  that  c  is  independent  of  x ,  and  thus  can 
be  dropped  from  the  test.  For  uncorrelated  noise,  our  model 
yields 


Grlrl(kfo>  =  Gss<kfo>  *  Gnlnl<k£o> 
Gr2r2<kfo>  =  +  Gn2n2<kfo> 

Grlr2<kfo>  =  “Gss<kfo>*-jk2'tf°I 


(2.11) 

(2.12) 

(2.13) 


(substituting  the  hypothesized  delay,  x,  for  the  true  delay, 
D).  Then  the  determinant  of  COV^j ^  in  (2.8)  is 
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o  ]  [aG  e 
J  1  ss 


+  jk2nf  t  n 
o  ] 


The  exponentials  in  the  last  term  cancel,  so  the  determinant 
and  thus  c  are  independent  of  x . 

Next,  we  observe  that  the  summand  of  (2.9)  is  ^0 
(positive  semidefinite,  [17,  p.  182]).  With  this  in  mind 
and  dropping  the  c  term,  the  equivalent  test  is  to  choose  x 
which  gives  the  minimum  value  for  in  (2.9). 

We  now  make  an  approximation  based  on  the  relationship 
between  the  Fourier  Series  and  the  Fourier  Transform.  For 
large  T  and  kfQ  held  constant,  it  can  be  shown  [18,  pp. 


23-25] 


=  SZ  R(f)*'  Q_1(f)  R(f)  df 


(2.14) 


We  next  determine  Q  ^(f)  from  (2.7),  substitute  into  (2.14), 
and  use  (2.11),  (2.12),  (2.13)  to  obtain 


J 


1 


J 


2 


J 


3 


where 


(f)  +  |r2( f )|  2/G 


(2.15) 


J3  =  2/“R1(f)R2(f)* 


df 


(2.16) 
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and 


'rlr2 


(f)  = 


Grlr2(f ^ 


Grlrl<f>Gr2r2<f> 

is  the  magnitude  squared  coherence  [19, 
bounded  by  -1  and  1. 


P-  54] 


(2.17) 
which  is 


Since  j£  is  not  dependent  on  t,  the  ML  estimate  reduces 
to  picking  t  which  maximizes  . 

If  the  R^(f)R2(f)  term  in  (2.16)  is  computed  for  the 
interval  (HT  by  averaging  short  segments,  then  it  can  be 
viewed  as  being  equal  to  T  times  an  estimate  of  the 
cross-power  spectrum,  T  Grlr2(f).  We  can  rewrite  J3  as 


J3  =  2T“  L 


Grlr2(f) 


G  (f) 
ss 


J2nft 


df 


;rlrl(f)Gr2r2(f)  ^Crlr2^ 


(2.18) 


but,  using  (2.11),  (2.12),  (2.13)  and  (2.17) 


Gss<f> 


;rlrl(f>Gr2r2(f> 


Grlr2(f ^ 


Grlrl<f>Gr2r2<f> 


Grlr2(f * 


Grlr2(f) 


Grlr2(f) 


Grlrl(f )Gr2r2(f ^ 
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Crlr2(f  ?. 


Grlr2<f> 


(2.19) 

Substituting  (2.19)  into  (2.18)  and  dropping  the  2T« 
constant,  we  get  our  ML  estimator  in  our  final  form,  that  is 
a  Generalized  Cross-Correlation  (GCC) 


}ML 

rlr2 


^  L»  Grlr2  ^ f 


Crlr2(f^ 


Grir2(f)|  [1  -  Crlr2(f)] 


-ei2nU  df 


(2.20) 


Summarizing,  our  ML  estimate  reduces  to  picking  as  our 
estimate  of  the  time  delay  D,  the  value  of  t  which  maximizes 
the  GCC.  The  GCC  is  computed  by  taking  the  inverse  Fourier 
Transform  of  a  weighted  estimate  of  the  cross-power  spectrum 
of  the  two  received  signals  r1(t)  and  r2(t).  The  ML 
weighting  is 

1  Crlr2<f> 

£)  - - 

|Grlr2<f)l  tl-Crir2<f>]  (2.21) 


which  is  defined  as  long  as  either  or  both  G  ,  n(f)  and 

nlnl 

Gn2n2^f^  are  non-zero  and  independent  (thus  C  ^ 2(f)  <1),  as 
would  be  the  case  in  practice. 


Up  to  this  point,  we  have  assumed  that  the  signal  and 
noise  auto-  and  cross-spectra  were  known.  In  practice  they 
may  have  to  be  estimated  by  any  of  various  techniques  (see 
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for  example  [19],  [20],  [21],  [22],  and  [23]).  In  such  case 

the  ML  estimator  is  only  approximately  achieved. 

We  now  will  very  briefly  show  that  this  ML  estimator 
has  a  variance  which  is  equal  to  the  minimum  bound  on  the 
variance,  and  as  such  is  optimum. 

Carter  and  Knapp  [7]  indicate  that  by  using  results 
obtained  from  [24],  the  variance  of  a  time  delay  estimate 
which  is  in  the  region  of  the  true  delay  (i.e.,  small 
estimation  error)  and  using  any  general  weighting,  W(f),  is 

£  w<f>2  (2nf)2Grlrl(f)Gr2r2(f)tl-Crlr2(f)]  df 

Var  [D]  =  - - - - - 

[T £  (2nf )  | Grlr2 ( f ) |  W( f )df ] 2 

(2.22) 

If  we  substitute  WML(f)  from  (2.21),  make  use  of  the 
definition  of  Crlr2(f)  from  (2.17),  and  consider  the 
frequency  symmetry  of  the  magnitude  of  the  cross-spectrum, 
(2.22)  reduces  to 


Var[D]  = 


2T  j“(2nf)‘ 


Crlr2<f> 


1-Crlr2^ 


df 


-1 


(2.23) 


We  now  show  that  this  is  the  minimum  variance.  The 
Cramer-Rao  Lower  Bound  (CRLB)  on  the  variance  is  [25,  p. 
72]  : 
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Var [D]  S  - 


a2ln  p(R(k)|  Q,  t ) 


ax' 


x  =  D 


Since,  for  the  ML  estimator,  the  only  part  that  is  dependent 
h< 

2 


on  x  is  J^,  then  we  can  write 


a 


Var [D]  £  - 


3x‘ 


E  ( J3/2 ) 


-1 


Since 


Grlr2<f>  *  Grlr2<f>  e 


-j2rtfD 


and 


E[Grir2(f)]  V*0  (f). 


rlr2 


we  have 


E  [J  /2]  =  T  j“e+j2nf(t-D)  . 


Crlr2<f> 


-  df 


1-Crlr2^^ 


Then  the  minimum  variance  is 


Min  Var[D]  = 


T/“(27tf)2 


Crlr2<f> 


1_Crlr2  ^  ^ 


df 


-1 


comparing  this  with  (2.23)  and  considering  the  symmetry  of 
Crlr2 ^ f ^  we  see  ^at  the  GCC  attains  this  minimum  variance. 
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2.3  Performance  of  the  GCC  Time  Delay  Estimator 

The  variance  of  the  TDE  by  the  GCC  method  (2.23),  as 
mentioned  earlier,  achieves  the  Cramer-Rao  Lower  Bound 
( CRLB )  of  the  estimator  variance.  In  the  form  of  (2.23),  it 
is  somewhat  difficult  to  interpret  directly  the  effects  of 
such  important  parameters  like  signal  and  noise  bandwidth 
and  SNR.  Quazi  [26]  gives  several  forms  which  are  more 
useful  for  such  interpretation. 

If  the  autospectra  of  the  signal  and  noise  are  assumed 
to  have  the  same  shape  (i~e.,  SNR  constant),  over  the  band 
fl  to  f2  HZ/  it:  is  easy  to  show  that  the  CRLB  is 

1 

f\  -  fl  (2.24) 


1 

4  -  (2.25) 

This  can  be  expressed  in  terms  of  signal  and  noise  bandwidth 
W  =  f2-f1,  and  center  frequency  f  =  (f  +f  )/2,  as  (for  SNR 
<<1) : 

1111  1 
_  *  _  •  •  • 

8zt2  SNR2  TW  fQ  l+(W2/12fQ3)  . 


1+2SNR 


CRLB  = 


8n2T  SNR2 


If  SNR  <<1  then 


CRLB  = 


8rt2T  SNR2 


CRLB 
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for  SNR  >>1 


3 


1 


1 


CRLB 


SNR 


or 


1 


111 


1 


CRLB 


From  the  above  representations,  it  is  seen  that  CRLB  varies 

inversely  as  the  time-bandwidth  product,  TW.  For  W<<f  ,  it 
,  o 

varies  inversely  as  the  center  frequency,  f  .  And  at  low 

SNR's,  it  varies  inversely  as  the  SNR2,  while  at  high  SNR's, 

it  varies  inversely  as  SNR.  Finally,  an  important 

relationship  is  that  for  a  given  SNR,  W,  and  fQ,  the 

variance  varies  inversely  as  the  total  observation  time  T. 

The  CRLB  gives  actual  achievable  performance  for  the 
GCC  if  erroneous  estimates  are  in  the  region  of  the  true 
correlation  peak.  If  there  is  insufficient  observation  time 
for  a  given  SNR  and  bandwidth,  then  large  anomalous  peaks 
can  occur  far  from  the  true  delay  region  which  results  in 
large  estimation  errors.  In  this  situation,  the  peformance 
of  the  estimator  will  be  signficantly  worse  than  the  CRLB. 
Figure  6  illustrates  this  effect  for  a  low  pass  filtered 
signal.  The  top  plot  shows  what  an  ideal  correlation  would 
look  like  for  no  noise  present  and  very  long  observation 
time.  For  the  realistic  case  of  moderate  noise  present,  the 
center  figure  is  more  representative.  The  correlation  peaks 
occur  away  from  the  true  delay,  but  are  still  in  the  general 
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ANOMALOUS  PEAKS  IN  THE  CROSS -CORRELATION  FUNCTION 
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neighborhood.  In  this  situation  the  CRLB  gives  a  good 
estimate  of  actual  performance,  i.e.,  it  is  a  "tight"  bound. 
As  the  SNR  is  decreased  further,  a  point  is  reached 
where  large  peaks  ( referred  to  as  anomalies  [27])  occur  far 
from  the  region  of  the  true  value  of  the  time  delay.  In 
this  situation,  the  performance  deviates  significantly  (gets 
worse)  from  the  CRLB,  and  it  is  no  longer  a  tight  bound. 


Recent  work  by  Ianniello,  Weinstein  and  Weiss  has  dealt 
with  determining  a  tighter  bound  than  the  CRLB  at  these 
lower  SNR's  [28],  [29],  [30],  [31].  Ianniello  [29]  has 

derived  a  performance  estimator  (not  a  bound  in  the  strict 
sense),  which  has  been  called  the  Correlation  Performance 
Estimate  (CPE)  [32],  for  the  GCC.  It  is  shown  to  be 

Var[°]CPE  =  PV3  +  U-p>vm[D)crlb 


is  the  correlation  peak  search  window,  i.e.,  if  we  know  a 
priori  that  the  maximum  possible  value  of  t  is  ±  Tq  then  we 
search  for  peaks  only  in  this  region  of  the 
cross-correlation  function  (for  example  if  we  know  the 
sensor  separation,  L,  we  know  the  maximum  delay  between  them 
is  Tq  =  ±  L/c,  c  =  speed  of  sound).  The  probability  that  an 
anomalous  noise  peak  occurs  (i.e.,  the  maximum  correlation 
peak  in  -Tq  ^  x  S  Tq  occurs  "far"  from  the  true  delay)  is  P. 
This  probability  is  shown  to  be: 


P  s  1-J 

J  r 


exp  [ -  ( x-K1 )  2/2  ]  dx  x  ( l/2n ) exp  ( -y2/2  ) dy^M’, 


2n 
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yBsT'(SNR) 

where:  =  - — 

[(SNR)2  +  (1+SNR)2]1/2 

SNR2 

Z  =  [1+ -  J1/2 

( 1+SNR)2 

M  =  4BTq  (for  uniform  spectra  over  -B^fgB) 

A 

=  #  independent  values  of  R  ,  _(t)  in  -T 

yiyz  o  o 

and  B  is  the  statistical  bandwidth  (B  =  2B  for  flat 

=>  S 

spectrum,  low-pass  signal  and  noise).  Uncorrelated 
zero-mean  Gaussian  signal  and  noise  are  assumed.  One  can 
evaluate  P  numerically  [13]. 

From  (2.26),  we  see  that  if  P  =  0  then  we  just  have  the 

CRLB.  If  P  =  1,  then  we  have  the  variance  of  a  uniformly 

distributed  random  variable,  ranging  between  ±  T  ,  as  might 

be  expected,  since  at  very  low  SNR's,  P  -»  1,  and  the  peak 

can  occur  anywnere,  with  equal  likelihood,  in  the  ±  T  peak 

o 

search  region. 

Scarbrough,  Tremblay  and  Carter  [32],  have  done 
computer  simulations  which  substantiate  these  results  and 
extend  them  to  make  some  observations  regarding  the 
significance  of  TDE  processing  by  coherent  vs.  incoherent 
means.  Figure  7  shows  a  typical  comparison  of  the  CRLB  and 
the  CPE  performance  plots.  (Note  that  the  vertical  axis  is 

A  A 

in  log1Q  of  the  standard  deviation  of  D,  Stdv[D],  not  the 
variance.  In  this  case  signal  and  noise  are  low-pass 
filtered,  flat  spectrum,  Gaussian  processes,  both  with 


FIGURE  7 
CPE  vs.  CRLB 
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bandwidth,  B  =  100  Hz,  integration  time,  T  =  8  secs,  and 
correlation  search  peak  windows  of  Tq  =  ±  1/8  sec.  The 
estimator  is  seen  to  achieve  the  CRLB  for  SNR  ^  6dB.  As  the 
SNR  is  lowered,  performance  rapidly  starts  to  deviate  away 
from  the  CRLB,  first  through  a  transition  region,  then 
finally  into  the  region  where  the  delay  estimate  becomes 
uniformly  random  across  the  peak  search  window,  ±  T  . 

Figure  8  is  a  family  of  curves  for  various  integration  times 
and  fixed  bandwidth  and  TQ.  As  the  integration  time  is 
increased  the  threshold  SNR  gets  lower.  Observe  that  for 
fixed  B  and  TQ,  the  value  of  log^Q  Stdv[D]  at  the  various 
threshold  SNR's  is  nearly  constant.  (Some  further 
discussion  on  this  is  given  in  [13]  and  [33] ).  Figure  9 
shows  the  results  of  the  computer  simulation  done  in  [32], 
for  two  integration  times  (the  X's  and  O's  are  the 
simulation  data  points  overlayed  on  the  calculated  CPE’s  and 
CRLB ' s ) . 

Figure  10  addresses  the  implications  regarding  coherent 
vs.  incoherent  TDE  processing.  The  top  curve,  labeled  "CPE 
(T=2)",  is  the  CPE  for  a  total  integration  time  of  2 
seconds.  If  4  of  these  delay  estimates  are  averaged  to  get 
a  new  estimate,  the  new  estimate  will  have  a  reduced 
variance  equal  to  1/4  that  achieved  by  2  second  estimations . 
This  process  is  referred  to  as  "incoherent"  processing  and 
its  expected  performance  is  the  curve  labeled  "CPE  (T=8, 

N=4)  (incoherent)."  If,  on  the  other  hand,  the  full  8 
seconds  of  data  was  processed  coherently  as  one  long 
correlation  (or  equivalently,  coherently  averaged 
consecutive  correlation  functions  from  contiguous, 
equal-sized  partitions  of  the  8  second  data),  the 
performance  obtained  is  the  curve  labeled  "CPE  (T=8)."  We 
see  that  at  higher  SNR's  the  8  second  incoherent  processing 


35 


'h 

o 

3 


SNR  (dB) 


FIGURE  8 

CPE  AND  CRLB  FOR  VARIOUS  INTEGRATION  TIMES 
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SNR  (dB) 


FIGURE  9 


CPE  AND  CRLB  VS.  SIMULATION  RESULTS 


FIGURE  10 

COHERENT  AND  INCOHERENT  INTEGRATION  PERFORMANCE 
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achieves  the  CRLB  for  T=8  seconds,  but  begins  to  deviate  at 
SNR2 ,  the  threshold  SNR  for  2  second  averaging.  But,  the  8 
second  coherent  processor  achieves  the  8  second  CRLB  down  to 
SNRg,  the  8  second  threshold  SNR  which  for  this  case  is  >4dB 
lower  than  SNR2.  Thus,  coherent  processing  is  optimum  at 
lower  SNR's  than  the  incoherent  method,  given  the  same 
amount  of  data. 


3 .  TIME  DELAY  ESTIMATION  WITH  MOVING 
SOURCE  AND/OR  RECEIVERS 


We  have  seen  that  the  TDE  correlator  performance  losses 
due  to  decreasing  signal-to-noise  ratio  (SNR)  can  be 
regained  to  near  CRLB  performance  by  increasing  the  coherent 
integration  time,  T.  For  source  and  receivers  with  fixed 
relative  positions,  the  GCC  can  be  used  without 
modification,  other  than  to  increase  the  total  integration 
time  (e.g.  by  taking  longer  segments  or  increasing  the  total 
number  of  segments  processed) .  However,  when  the  relative 
motion  between  the  source  and  receivers  is  different  for 
each  of  the  receivers,  the  time-delay  of  arrival  varies  over 
T.  This  results  in  a  loss  of  coherence  between  the  2 
received  signals  which  worsens  as  T  gets  longer,  effecting  a 
"smearing"  of  the  correlation  peak  and  a  drop  in  it's 
amplitude  which  in  turn  can  cause  a  significant  increase  in 
error.  This  is  especially  true  if  the  delay  varies  by  more 
than  the  correlation  time  width  of  the  source  signal  during 
the  [0+T]  observation  time  [34] .  In  order  to  get  improved 
performance,  one  must  compensate  for  this  varying  time  delay 
to  again  allow  longer  integration  time  and  use  of  the  GCC 
processor.  This  section  addresses  this  problem  and  some 
approaches  to  resolving  it. 

3.1  Motion  Model 


The  model  presented  here  is  based  on  the  one  given  by 
Knapp  and  Carter  [35].  Other  versions,  which  are  eqivalent 
but  whose  form  offer  additional  insight  to  the  problem  are 
presented  in  [34],  [36],  [13],  and  [38].  The  form  of  the 

model  in  [13]  is  primarily  used  in  this  section. 
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We  let  s'(t)  be  the  signal  received  at  receiver  1  when 
no  relative  motion  or  noise  is  present.  Then,  for  the  case 
with  relative  motion  and  additive  noise,  an  appropriate 
model  would  be 

rx(t)  =  s’  (p^tj+n^t)  (3.1. a) 

r2(t)  =  «s'(p2(t+D))+n2(t)  (3.1.b) 

The  functions  r^t)  and  r2(t)  are  the  receiver  1  and  2 
inputs  having  stationary,  white,  zero-mean  additive  noise 
n1(t)  and  n2(t),  which  are  uncorrelated  with  each  other  and 
s'(t).  The  time  delay,  D,  is  the  time  difference  of  arrival 
(TDOA)  as  in  the  no  motion  model,  and  a  is  the  attenuation 
of  the  signal  at  receiver  2  relative  to  that  at  receiver  1 
(in  most  cases  of  interest,  asi).  The  relative  motion 
between  the  source  and  each  of  the  receivers  causes  a  time 
compression/expansion,  termed  "time  companding,"  [35]  which, 
in  general,  is  different  for  each  receiver,  and  is 
represented  by  p^  and  p2.  For  relative  velocity  v  between 
the  source  and  receiver  i,  and  signal  propagation  velocity, 
c ,  then 


P-L  =  l+v1/c 

P2  =  1+v2/c 

For  acoustic  sources  in  the  ocean,  c  s  4900  feet/second  and 
generally  c>>v^,  so  p^  si.  If  a  positive  velocity  v^  is 
assigned  to  source  and  receiver  motion  towards  each  other, 
then  p^  >  1.  It  can  be  seen  that  the  model  in  (3.1)  reduces 
to  the  no  motion  model  for  P1  =  P2  =  1.  An  important 
assumption  made  is  that  p^  and  p2  are  essentially  constant 
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over  the  correlation  time  0->T  (i.e.,  accelerations  in 
source/receiver  positions  are  ignored). 

A  slightly  different  representation  of  the  model  will 
give  a  little  added  insight  to  the  problem. 

Let  s ( t )  =  s*(P1t).  Then 

s' (t)  =  s(t/p1) 

s'(P2(t+r)))  =  s([p2/p1](t+D)) 


Now  let 


e  =  i2/»x 

then  we  can  write  (3.1)  as 

rx ( t)  =  s(t)  +  nx  (t)  (3. 2. a) 

r2(t)  =  ocS(p(t+D))  +  n2(t)  (3.2.b) 

In  this  form  we  see  that  the  motion  effects  can  be  thought 
of  in  terms  of  a  relative  time  companding  ( RTC ) ,  [35],  [38], 

of  the  signal  at  receiver  2  with  respect  to  that  at  receiver 

1. 


In  yet  another  form,  the  motion  effects  can  be  seen  as 
a  time-varying  time  delay.  From  (3.2.b): 
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s(p(t+D))  =  s(pt  +  pD) 

=  s(t-t+pt+pD) 

=  s(t+(p-l)t+pD)  (3.3) 


Define : 


d'  =  (P-1) 
dQ  ^  PD 

d(t)  ^  d't  +  do  (3.4) 


then 


s ( p ( t+D ) )  =  s(t+d't+d  ) 

=  s ( t+d ( t ) ) 

and  our  model  can  now  be  written  as : 

rx(t)  =  s(t)  +  nx(t)  (3. 5. a) 

r2(t)  =  «s(t  +  d( t) )  +  n2(t)  (3.5.b) 

where  d(t)  is  seen  to  be  a  linearly  time-varying  delay  with 
rate,  d’,  and  initial  value,  dQ.  The  two  representations  of 
the  motion  model  in  (3.2)  and  (3.5)  point  to  two  methods  for 
compensating  for  motion  effects.  The  first  indicates  a  need 
to  compensate  by  processing  the  receiver  2  inputs  with  a 
time  compander  that  negates  the  time-companding  due  to 
motion  and  reverts  back  to  the  no  motion  case.  The  second 
form  indicates  compensation  by  removal  of  the  time-varying 
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time  delay  effects.  These  approaches  are  discussed  in  the 
next  two  sections . 

3 . 2  Compensating  for  Relative  Time  Companding 

From  (3.2)  it  is  apparent  that  the  effect  of  relative 
time  companding  may  be  removed  by  passing  r2(t)  through  a 
time-compander  which  will  perform  time-scaling  that  is  the 
inverse  of  p,  i.e.  generate 

r3(t)  =  r2(t/b)  (3.6) 

where  b  =  p,  ideally.  Then  r^(t)  and  r3(t)  can  be  processed 
with  the  GCC  as  in  the  no-motion  case. 

By  an  argument  very  similar  to  that  used  to  derive  the 
GCC  as  a  ML  estimator,  Knapp  and  Carter  [35]  show  that  the 
ML  estimator  when  motion  is  present  (at  least  to  a 
reasonable  approximation  when  the  power  spectra  are 
unknown),  is  to  maximize 

J(t/b)  =  R1(f)R3(f,b)*W(f,b)e:i27tftb  df  (3.7) 

where  x  and  b  are  hypothesized  values  for  D  and  p 
respectively,  R-^(f)  and  R3  (f,b)  are  the  Fourier  Transforms 
of  r1(t)  and  r3(t),  and  W(f,b)  is  the  ML  weighting  function 
derived  from  estimates  of  the  auto-  and  cross-power  spectra 
of  r1(t)  and  r3(t),  defined  as 


W(f,b) 


crir3(f'b> 


Grlr3(f'b> 


( i-c 


rlr3 


(f,b))  (3.8) 
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and  Crir3(f/b)  is  the  magnitude  squared  coherence  estimate: 


C(f,b)  = 


lGrlr3(f'b> 


Grlrl<f>  Gr3r3<f'b> 


(3.9) 


Note  that  except  for  the  dependencies  on  b  due  to  the 
hypothesized  compensation  of  ^(t)  the  above  are  the  same 
form  as  for  the  GCC.  Then  J(x,b)  can  be  seen  to  be  a 
correlation  function  with  argument  xb,  which  if  t  =  D  and  b 
=  0  is  Dp  =  dQ. 


The  compensation  for  time-companding  must  be  done  for  a 
range  of  hypothesized  b's,  then  each  J(x,b)  is  computed  as 
in  (3.7)  giving  a  two-dimensional  ambiguity  surface  in  t  and 
b  space.  The  global  peak  on  this  surface  is  found  and  it's 
corresponding  values  of  x  and  p  are  our  ML  estimates  of  D 
and  p  respectively.  Figure  11  illustrates  this  process.  The 
COMP^,  i  =  l  to  Nc,  represents  compensation  for 
time-companding  using  one  of  Nc  hypothesized  values  for  b. 

The  preceding  process  is  quite  expensive 
computationally .  The  entire  GCC  process  must  be  done  for 
each  hypothesized  b,  including  calculation  of  W(f,b).  The 
compensating  companding  processor  could  involve  hardware 
such  as  a  variable  speed  analog  tape  recorder  or  multiple 
rate  samplers  to  effect  the  time  compression/expansion  of 
r2 ( f ) •  0r  it  could  be  done  in  software  with  a  very  fine 
resolution  interpolation  algorithm  followed  by  variable  rate 
re-sampling.  Although  any  of  these  methods  are  feasible 
(Betz  has  implemented  the  latter  technique  and  found  it  to 
perform  near  optimum  [39]),  it  would  be  desireable  to  find  a 
more  efficient  way  to  compensate  for  the  effects  of  relative 
motion.  As  will  be  seen  in  the  next  section,  this  can  be 
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FIGURE  11 

TIME -COMPANDING  BLOCK  DIAGRAM 
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done  if  we  compensate  for  the  varying  time-delay  of  (3.5) 
instead  of  the  relative  time  companding  of  (3.2). 

3.3  Compensating  for  Varying  Time-Delay 

The  general  strategy  that  will  be  used  to  compensate 
for  the  linearly  time-varying  delay  indicated  in  (3.4)  and 
(3.5)  can  be  summarized  briefly  as  follows.  The  total 
coherent  processing  interval  from  t  =  0  to  T  seconds,  is 
divided  into  several  equal-sized,  short-time  intervals. 
Assuming  that  the  delay  change  during  each  short-time 
interval  is  relatively  small,  a  modified  model  is  formed  by 
approximating  the  continuously  changing  delay  during  each 
short-time  interval  as  a  constant  equal  to  the  actual  delay 
at  the  mid-point  of  that  interval.  Then,  based  on  this 
model,  the  compensation  is  done  in  a  step-wise  fashion  for 
each  short-time  interval,  by  adjusting  for  the  amount  of 
delay  change  from  interval  to  interval.  A  compensated 
short-time  GCC  type  cross-correlation  function  is  computed 
for  each  interval.  These  are  then  averaged  to  get  a  total 
coherent  processing  time  of  T  seconds.  The  above  is  done 
for  a  range  of  assumed  delay  rates,  forming  an  ambiguity 
surface  from  which  the  ML  delay  and  delay  rate  are 
determined.  This  process  closely  approximates  the  ML  method 
of  relative  time-companding  compensation  of  Section  3.2,  but 
without  the  need  of  performing  time-scaling.  Also,  only  one 
weighting  function,  W(f),  needs  to  be  computed  since  there 
is  no  dependency  on  an  assumed  time-scaling  factor,  b,  as 
was  the  case  in  the  time-companding  compensation  method. 

This  piecemeal  approach  to  removing  the  effects  of 
varying  time-delay  can  be  implemented  in  two  ways,  as 
described  in  [36]  and,  more  recently,  by  [13],  [37]  and 
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[38].  The  effect  of  this  process  has  been  described  as 
"deskewing"  in  [38]  for  reasons  that  will  be  apparent 
shortly. 

From  (3.5)  we  can  write  the  cross-correlation  function 
of  r1(t)  and  r2(t)  at  times  and  t2  as: 

Rrlr2  =  Etr1(t1)  r2(t2)] 

=  e{  [sftj)  +  n1<t:1 )  ]  [«s(t2  +  d(t2))  +  n2(t2)]J 

Since  s(t),  n^(t)  and  n2(t)  are  zero-mean  and  mutually 
uncorrelated,  we  have 

Rrlr2^tl't2^  =  ccE[s(t1)s(t2  +  d(  tp  ) )  ] . 

Letting  =  t  +  x  and  t2  =  t, 

Rrlr2(t+T,t)  =  aE[s(t+T )s(t+d(t) )] 

=  «Rss(t+x,t+d(t)) 

Since  s(t)  is  wide  sense  stationary,  its  auto-correlation 
function  is  dependent  only  on  the  difference  of  its  time 
arguments ,  so 

Rrlr2(t+T't)  =  aRSS(T-d(t)> 

=  aRss(t-d't-d0)  (3.10) 

We  see  that  the  cross-correlation  function  is  dependent  on 
t,  so,  clearly,  r^(t)  and  r2(t)  are  jointly  non-stationary . 
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Our  next  step  will  be  to  break  up  the  desired  coherent 
processing  time,  t  =  0  T,  into  M  equal-sized  short  time 
intervals  and  label  them  Tk  where 


Tk  = 


t  ( k-1 )T/M  ^  t  ^  kT/M 


k  =  1,2,  ...  ,M 


It  is  desired  to  make  these  intervals  short  enough  so  that 
the  amount  of  change  in  the  delay,  d(t),  is  small  compared 
to  the  change  over  the  entire  interval.  Yet,  they  should  be 
long  enough  to  assure  independent  correlations  for  each 
interval.  Betz  [38]  has  numerically  determined  an  upper 
bound  on  the  number  of  segments  to  be 


M  ^  WT/3.5 


for  low-pass,  white  signal  and  noise. 

Next,  we  approximate  the  continuous,  linearly  varying 
time  delay,  d(t),  in  a  piece-wise  constant  fashion  by  using 
its  values  at  the  midpoint  of  each  interval,  T^,  [13],  [36], 

[38] .  The  value  of  t  at  such  midpoints  is  denoted 


tk  =  ( k-1/2 )T/M 

and  the  value  of  the  delay  at  the  midpoint  is 


d(tk)  =  d'tk  +  dQ 


=  dk  +  do 


(3.11) 


where  d. 


=  d't. 


k* 


This  is  illustrated  in  Figure  12 . 
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STEP-WISE  CONSTANT  APPROXIMATION  OF  DELAY  RATE 
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Substituting  (3.11)  into  (3.10)  gives  the  approximate 


time  cross— correlation  function  for  each  interval 
as 


T 


k 


Rrlrl(T'k>  =  *Rss(T-  dk  -  do>  (3  •  !2 ) 

which^ow  is  wide  sense  stationary  (independent  of  t)  within 
the  k  interval.  In  this  form,  we  can  envision  that  over 
the  total  observation  time,  t  =  0  ->  T,  that  the  short  term 
cross-coi relation  functions  are  approximately  equal  to 
short-time  auto-correlation  functions  of  the  signal  which  is 
initially  shifted  in  delay  by  dQ,  and  subsequently  shifted 
(or  "skewed")  in  x  by  the  amount  dk  during  each  following 
short-time  interval.  Figure  13  is  an  actual  set  of  these 
short-time  cross-correlations  generated  by  the  simulation 
program  described  in  Appendix  A,  which  illustrates  this 
skewing  effect  for  dQ  =  0,d’ =1/256,  and  SNR  =  +6dB.  The 
plots  show  an  expanded  region  around  +50  sample  intervals . 

It  should  be  clear  from  the  figure  that  if  these 

lutions  are  simply  averaged,  without  compensating  for 
the  skewness,  that  the  resulting  cross-correlation  function 
will  be  smeared,  as  shown  in  Figure  14.  This  smearing  leads 
to  a  broadening  of  the  correlation  peak  and  a  reduction  in 
it's  mean  value  [34],  [36],  [38].  This  can  be  expected  to 
increase  the  uncertainty  of  the  true  peak  location  and 
susceptibnty  to  random  noise  effects,  resulting  in  an 
increase  in  the  variance  of  the  delay  estimate.  Another 
effect  is  to  cause  a  bias  in  the  delay  estimate  towards  the 
value  of  the  true  delay  at  time  T/2  [35].  So  our  intent  is 
to  compensate  by  realigning  or  "deskewing"  [38]  the 
correlation  functions  with  respect  to  that  of  the  beginning 
of  the  T^  interval  before  summing,  in  order  to  get  the 
maximum  benefit  of  the  averaging  process  in  reducing  noise. 
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FIGURE  13 

UNCOMPENSATED  SHORT-TIME  CROSS -CORRELATION 
FUNCTIONS  (UNAVERAGED) 
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FIGURE  14 

UNCOMPENSATED  SHORT-TIME  CROSS-CORRELATION 
FUNCTIONS  (AVERAGED) 
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As  will  be  seen  next,  this  deskewing  can  be  done  as  a  phase 
rotation  of  the  short-time  cross-spectra  in  the  frequency 
domain,  or  as  a  delay  shift  in  the  time-delay  domain.  The 
first  method  has  been  implemented  by  Kuhn,  Robison,  and 
Winfield  [36],  by  Scarbrough  [13],  and  is  the  method  used  in 
the  simulations  done  for  this  thesis,  described  in  Section 
4.0  and  Appendix  A.  The  latter  method  has  been  implemented 
by  Betz  [38],  and  termed  "Deskewed  Short-Time  Correlation," 
an  name  which  could  be  appropriate  for  either  approach,  as 
the  net  effect  is  essentially  the  same. 

3.3.1  Compensation  in  the  Frequency  Domain 
(Phase  Rotation) 

We  proceed  by  taking  the  Fourier  Transform  of  (3.12)  to 

get 

Grlr2(f'k)  =  “Gss(f)e"j27tf(dk+do) 

It  is  clear  that  the  time-varying  delay  is  equivalent  to  a 
corresponding  time-varying  phase  rotation  of  the  signal 
auto-power  spectrum  by  an  amount  (-2nfdk)  radians.  We  can 
compensate  for  this  with  a  counter  rotation,  +  2nfdk,  where 


is  the  amount  of  delay  shift  at  time  t,  (since  t  =  0);  t '  is 
an  hypothesized  delay  rate  (and  so  an  estimate  of  the  true 
rate,  d * ) ,  chosen  from  a  range  of  possible  delay  rates. 

Then,  averaging  the  compensated  short-time  cross-spectra: 
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M 


Grir2  (f'T)  =  (1/M)  Grlr2(f'k)e+j27t(T'tk)  (3'13) 


M 


=  ( l/M )  £  a  Gss(f,k)e-:i27Tf(dk-dk  +  do} 


k-1 


If  we  hypothesize  x',  and  thus,  d^,  correctly  then 
Grlrl  <f>  =  «Css(f)e-i2Rfio 

This  is  the  same  cross-spectrum  used  in  the  GCC  estimator 
(2.20),  so  that  after  delay  rate  compensation  the  no-motion 
GCC  processing  (i.e.  ML  weighting  and  inverse  Fourier 
Transform)  can  be  used  to  estimate  the  delay  for  a  range  of 
hypothesized  t'.  The  resultant  two-dimensional  (x,  x') 
ambiguity  surface  is  searched  for  the  global  maximum  peak. 
The  corresponding  value  of  x  and  x '  are  our  ML  estimates  for 
d'  and  dQ  respectively.  This  process  is  illustrated  in 
Figure  15. 


3.3.2  Compensation  in  the  Delay  (Deskewing) 

In  this  approach,  the  individual  cross-spectra  from 
each  short-time  interval  is  weighted  by  W^^(f)  and  inverse 
Fourier  Transformed  to  get  M  short-time  cross-correlation 
functions.  Each  of  these  is  progressively  skewed  in  the 
delay  dimension,  as  in  Figure  13.  Then  each  short-time 
correlation  function  for  the  interval  T^  is  shifted 
(deskewed)  by  d^  =  x't^.  This  is  done,  for  example,  by 
interpolation  with  a  cubic  spline  as  in  [38].  Finally  the 
deskewed  cross-correlations  are  averaged  to  get  the 
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FIGURE  15 

PHASE  ROTATION  COMPENSATION  BLOCK  DIAGRAM 
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equivalent  no-motion  GCC.  This  is  done  for  a  range  of 
hypothesized  t ' ,  again  generating  a  two  dimensional 
ambiguity  surface  from  which  the  estimates  for  d'  and  d 

o 

obtained.  This  process  is  illustrated  in  Figure  16. 

The  simulations  reported  in  [38]  verify  that  this 
method  very  nearly  approximates  the  ML  method  of  section 
3.2,  for  a  wide  range  of  d' ,  and  requires  an  order  of 
magnitude  fewer  computations  than  the  ML  method. 


are 
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FIGURE  16 

DESKEWING  COMPENSATION  BLOCK  DIAGRAM 


4.0  SIMULATION  OF  MOTION  COMPENSATED  TDE 


A  computer  simulation  of  the  motion  compensated 
cross-correlator  described  in  Section  3.3.1  was  implemented 
on  a  Hewlett  Packard  HP9845B  desktop  computer,  to 
experimentally  observe  the  effects  on  performance  when 
motion  compensation  is  used  or  not  used,  when  there  is 
relative  motion  present.  It  was  also  done  to  verify  some  of 
the  key  theoretical  analysis  presented  earlier,  and  to  gain 
some  insight  into  the  problems  associated  with 
implementation.  The  simulation  process,  test  run 
parameters,  and  test  results  and  discussion  are  covered  in 
the  following  subsections.  A  description  of  the  simulation 
program  is  provided  in  Appendix  A. 

4 . 1  Simulation  Description 

4.1.1  Signal  Generation 

Simulated  samples  of  the  receiver  signals  r^(t)  and 
r2(t)  °f  (3.5)  were  generated  as  follows.  Three  independent 
Gaussian  pseudo-random  noise  sequences,  having  zero-mean  and 
unit  variance,  were  generated  from  independent  uniform, 

[0-1] f  sequences  by  the  "direct  method"  of  [40,  pp.  953]. 

One  of  these  was  used  as  the  sample  sequence  s(n)  from  s(t). 
The  other  two,  n-^(n)  and  n2(n),  were  scaled  according  to  the 
desired  SNR,  and  represented  samples  from  n^(t)  and  n_(t). 

The  s(n),  n^n),  and  n2(n)  sequences  are  next  filtered 
with  a  Low-Pass  4-stage,  2-pole  Butterworth  HR  digital 
filter  with  fmax  =  100  Hz,  the  upper  frequency  of  the  pass 
band  for  an  assumed  sample  rate  of  2048  Hz. 
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The  s(n)  sequence  was  processed  by  a  time-varying 
filter  function  which  introduced  a  phase  shift  that  was 
linearly  dependent  on  frequency,  and  was  varied  from  time 
sample  to  time  sample,  according  to  the  desired  simulated 
delay  rate,  d'  (41).  The  frequency  dependent,  linear  phase 
shift  in  the  frequency  domain  is  equivalent  to  a  delay  in 
the  time  domain  [19,  pp.  66-67].  Varying  this  phase  for 
each  sample  simulates  a  time-varying  time  delay,  d(n),  that 
varies  from  sample  to  sample  at  the  rate  of  d'  sample 
intervals/sample  interval.  A  relative  attenuation,  «=1,  is 
simulated.  Finally  s(n)  and  s(n+d(n))  are  added  to  n^(n) 
and  ^(n)  to  get  r^(n)  and  ^(n),  the  simulated  samples  from 
r1(t)  and  ^(t).  Sequence  lengths  of  256  samples  were  used 
for  each  Tk  interval. 

4.1.2  Processing 

The  time-varying  time-delay  compensation  technique  of 
Section  3.3.1  was  implemented  as  follows.  Since  the  signal 
and  noise  generated  have  flat  spectra  and  equal  SNR  across 
the  band,  an  ML  weight  of  W^*(f)  -  \  was  used  (i.e.,  the  SCC 
weight) .  The  frequency  domain  cross-correlation  algorithm 
of  [21,  pp.  555-562]  was  used,  with  512  point  FFT's 
(including  256  augmenting  zeroes),  but  prior  to  inverse 
transforming,  the  cross-spectra  were  phase  rotated  according 
to  (3.13).  Five  hypothesized  values  for  x'  were  used: 
-2/256,  -1/256,  0,  1/256,  and  2/256  seconds/second  (or 
equivalently,  sample  intervals/sample  interval).  Each 
phase-rotated  cross-spectrum  was  then  summed  and  averaged 
with  the  corresponding  ones  from  the  previous  intervals  and 
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then  inverse  FFT'd  after  all  M  intervals  had  been  processed. 
The  five  resulting  cross-correlations  thus  formed  a  5  by  512 
point  ambiguity  surface  in  x ' ,  t  space.  The  global  peak  was 
determined  and  it's  corresponding  values  for  x  and  x'  were 
used  as  the  ML  estimates  of  dQ  and  d' ,  respectively. 

4.1.3  Simulation  Parameters 

The  simulation  was  conducted  using  a  total  coherent 
integration  time  of  T  =  2  seconds,  and  an  assumed  sample 
rate  of  2048  H  .  The  total  time,  T,  was  sub-divided  into 
sixteen  1/8  second  intervals  (i.e.  M=16)  of  256  points  each 
and  processed  as  indicated  above.  SNR’s  ranging  from  -12dB 
to  +6dB  in  3dB  steps  were  used,  with  100  trials  (2  seconds 
each)  run  at  each  SNR.  Delay  rates  of  d'  =  ±1/256 
seconds/second  and  initial  delays  of  dQ  =  +16,  -16  and  0 
samples  were  simulated.  Statistics  on  the  resulting  delay, 
delay  rate,  and  correlation  peak  amplitude  were  gathered  for 
both  compensated  and  uncompensated  correlations. 

4.1.4  Results 


Figure  17  shows  the  compensated  (by  phase-rotation) 
short-time  cross-correlations  for  each  short  segment  in  a 
single  trial.  The  SNR  was  +6dB  (this  is  the  same  trial 
illustrated  in  Figures  13  and  14). 

Figure  18  is  the  ambiguity  "surface"  generated  for  the 
same  trial  shown  in  Figures  13,  14,  and  17. 

Table  2  lists  the  error  performance,  log1Q  Stdv[d  -d  ] 
for  the  delay  estimates  for  compensated  and  uncompensated 
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FIGURE  17 

COMPENSATED  SHORT-TIME  CROSS-CORRELATIONS 
(PHASE  ROTATION  METHOD) 
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FIGURE  18 

SAMPLE  AMBIGUITY  SURFACE 
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SNR(dB) 

L°g10 

Stdv[d  -d  1 
o  o 

Mean  Peak  Loss 

COMP 

NO  COMP 

6 

0 

-3.41 

-1.63 

3 

-4.13 

-3.36 

-1.62 

0 

CD 

» 

CO 

1 

-3.28 

-1.60 

-3 

-3.46 

-2.52 

-1.58 

-6 

-1.59 

-1.39 

-2.25 

-9 

-1.22 

-1.21 

-2.72 

-12 

-1.16 

-1.17 

-4.16 

TABLE  2 

SIMULATED  RUN  STATISTICS 
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processing.  The  mean  peak  loss  is  defined  as  10  log1Q  of 
the  ratio  of  the  mean  correlator  peak  value  for  the 
uncompensated  processing  to  that  of  the  compensated 
processing. 

Figure  19  shows  plots  of  the  error  performance  for  the 
simulation  runs  compared  to  those  predicted  by  the  CPE  and 
CRLB.  The  x's  are  from  a  separate  simulation  reported  in 
[32],  for  which  there  was  no  relative  motion.  The  points 
marked  with  "dots"  are  from  the  simulation  with  motion 
reported  here.  Note  that  the  vertical  axis  is  in 
l°9l0  Stdv[dQ-d0].  The  confidence  intevals,  [20,  pp. 
113-115],  are  indicated.  The  upper  plot  shows  the 
uncompensated  performance,  and  the  compensated  results  are 
shown  in  the  lower  plot. 

4.1.5  Discussion 


Figure  17  clearly  shows  the  effectiveness  in  the 
compensation  scheme  for  d' =1/256.  The  bottom  correlation  is 
the  average  of  all  16  segments.  The  peak  has  been  shifted 
back  to  the  true  value  of  dQ=0  as  expected  from  (3.13). 

Figure  18  shows  the  final  corss-correlation  from  a 
single  trial  for  each  hypothesized  t 'which,  together,  form 
the  ambiguity  surface.  The  center  one  corresponds  to  t'=0, 
or  no  comenpsation.  Note  the  difference  in  its  peak  width 
and  height  compared  to  the  correctly  compensated  correlation 
( t 1 =1/256),  second  from  the  bottom.  As  predicted  earlier, 
the  uncompensated  correlation  peak  is  both  broader 
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FIGURE  19 

COMPENSATED  VS.  UNCOMPENSATED 
PERFORMANCE ,  SIMULATION  RESULTS 
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("smeared")  and  lower  in  amplitude,  due  to  uncompensated 
averaging. 

The  performance  plots  of  Figure  19  are  consistent  with 
the  theory  presented.  The  threshold  SNR  for  the 
uncompensated  estimates  appears  to  occur  at  a  higher  SNR  (= 

+  1.5dB  higher)  than  for  the  compensated  correlations,  which 
is  consistent  with  the  s  1.6dB  loss  in  peak  value.  .Also, 
the  error  is  greater  at  all  SNRs  above  the  threshold  due  to 
spreading  and  drop  in  amplitude  of  the  uncompensated 
correlation  peak  which  is  expected  to  cause  an  increase  in 
local  error  variation. 


5.0  CONCLUSIONS  AND  RECOMMENDATIONS 


The  theory  and  analysis  presented  earlier  had  indicated 
that  the  GCC  is  a  near  optimum  (when  spectra  must  be 
estimated)  estimator  of  time  delay  for  fixed  source  and 
receivers.  It  was  shown  that  the  introduciton  of  motion 
would  require  compensation  for  the  effects  of  relative 
time-companding,  or  equivalently,  time  varying  delay; 
otherwise  a  reduction  in  performance  would  be  encountered. 
Several  schemes  for  compensating  for  these  effects  were 
described  and  one  using  phase  rotation  of  the  estimated 
cross-sectra  was  implemented  in  the  simulation  program 
described  in  Appendix  A.  The  results  of  running  100  trials 
for  a  range  of  SNR's  are  consistent  with  the  theory  and 
analysis  presented  earlier.  In  particular,  compensated 
estimations  performed  near  the  optimum  defined  for  the  case 
with  no-motion.  The  effect  of  no  compensation  was  an 
overall  increase  in  estimation  error  variance  and  an  onset 
of  anomalous  errors  at  a  higher  SNR  than  for  the  compensated 
case. 


The  variety  of  the  trials  run  for  the  simulation 
results  presented  here  was  very  limited  in  scope  due  to  the 
enormous  computation  times  involved.  Use  of  an  array 
processor  (such  as  reported  in  [32]),  would  greatly  expand 
the  number  possible  run  parameters  which  could  be  simulated. 

It  is  recommended  that  more  simulations  be  run  using  a 
wider  range  of  delay  rates.  Also,  simulation  results  are 
needed  which  use  a  larger  number  of  hypothosized  delay-rate 
compensations,  to  determine  the  effects  that  more  points  on 
the  ambiguity  surface  might  have  on  the  number  of 
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occurrences  of  anomalous  peaks.  Finer  SNR  steps  would  help 
to  better  define  the  shape  of  the  performance  curves, 
allowing  a  more  accurate  empirical  comparison  of  error 
performances . 
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APPENDIX  A 


Figure  A.l  is  a  program  structure  chart  indicating  the 
subroutines  of  the  simulation  program  and  their  relationship 
to  each  other.  Figure  A. 2  is  a  functional  flow  chart 
indicating  the  top-level  flow  of  the  simulation  program 
process . 

For  listings  or  program  copies  contact: 


R.  J.  Tremblay,  Code  3314 
Naval  Underwater  Systems  Center 
New  London,  Connecticut  06320 

or 

Professor  D.  W.  Lytle 
EE  Building,  FT-10 
University  of  Washington 
Seattle,  Washington  98195 


FIGURE  A. 1 

PROGRAM  STRUCTURE  CHART 
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